function [res,resEuler,Evar_cu] = ResidualEquationProjExpectedControl(thetaVec,setup)
 
% This function computes the residuals for the problem for computing E_t[var_cup] 
% We use the notation that Evar_cu denotes E_t[var_cup], meaning that the
% problem is -Evar_cu + var_cup = 0
 
%% Getting some information from setup 
scaleX_in_g = setup.scaleX_in_g;
out         = setup.out;
x_cu        = setup.xGrid;
numX        = size(x_cu,1);
nNodes      = setup.n_nodes;
wNodes      = transpose(setup.weight_nodes);
epsNodes    = setup.epsi_nodes;
params      = setup.params;
appOrder    = setup.appOrder;
nx          = setup.nx;
ny          = setup.ny;
nx1         = setup.nx1;
gSS         = setup.gSS;
hSS         = setup.hSS;
% The unknowns: for Evar_cu = E_t[var_cup]
thetaProj   = reshape(thetaVec,1,setup.nDim2Theta);
g0          = thetaProj(1,1); 
gx          = thetaProj(1,2:1+nx); 
if appOrder > 1 
    gxxV    = thetaProj(1,nx+2:1+nx+nx*(nx+1)/2); 
end 
if appOrder > 2 
    gxxxV   = thetaProj(1,1+nx+nx*(nx+1)/2+1:1+nx+nx*(nx+1)/2+nx*(nx+1)*(nx+2)/6); 
end 
% The knowns: for var_cu  
thetaVar    = reshape(setup.thetaVar,1,setup.nDim2Theta);
var0        = thetaVar(1,1); 
varx        = thetaVar(1,2:1+nx);
if appOrder > 1 
    varxx   = thetaVar(1,nx+2:1+nx+nx*(nx+1)/2);
end 
if appOrder > 2 
    varxxx  = thetaVar(1,1+nx+nx*(nx+1)/2+1:1+nx+nx*(nx+1)/2+nx*(nx+1)*(nx+2)/6); 
end 
 
%Unfold params 
BETTA= params.BETTA;
B= params.B;
CHI= params.CHI;
CHI0= params.CHI0;
THETA= params.THETA;
DELTA= params.DELTA;
ALFA= params.ALFA;
PHI= params.PHI;
PHIzero= params.PHIzero;
ZI= params.ZI;
KAPAw= params.KAPAw;
ETA= params.ETA;
RHOR= params.RHOR;
PHIpai= params.PHIpai;
PHIpai_1= params.PHIpai_1;
PHIy= params.PHIy;
PHIy_1= params.PHIy_1;
PHIc= params.PHIc;
PHIc_1= params.PHIc_1;
PHIl= params.PHIl;
NU= params.NU;
U0= params.U0;
U0d= params.U0d;
OMEGAz= params.OMEGAz;
OMEGAd= params.OMEGAd;
OMEGAn= params.OMEGAn;
OMEGAp= params.OMEGAp;
OMEGAa= params.OMEGAa;
RHOz= params.RHOz;
RHOd= params.RHOd;
RHOn= params.RHOn;
RHOp= params.RHOp;
RHOa= params.RHOa;
Kss= params.Kss;
OUTPUTss= params.OUTPUTss;
PAIss= params.PAIss;
MUZss= params.MUZss;
Css= params.Css;
AA= params.AA;
lss= params.lss;
Rss= params.Rss;
Wss= params.Wss;
 
%% Unfold the states and controls for the core of the model - all in levels 
% Getting the states at time t 
c_ba1 = out.c_ba1;
muz_cu = out.muz_cu;
d_cu = out.d_cu;
n_cu = out.n_cu;
paistar_cu = out.paistar_cu;
a_cu = out.a_cu;
 
% Getting the controls at time t 
c_cu = out.c_cu;
pai_cu = out.pai_cu;
evf_cu = out.evf_cu;
 
% Getting the states at time t+1 
c_ba1p = out.c_ba1p;
muz_cup = out.muz_cup;
d_cup = out.d_cup;
n_cup = out.n_cup;
paistar_cup = out.paistar_cup;
a_cup = out.a_cup;
 
% Getting the controls at time t+1 
c_cup = out.c_cup;
pai_cup = out.pai_cup;
evf_cup = out.evf_cup;
 
% Computing the new control variable  
x1_cu  = x_cu(:,1)/scaleX_in_g(1,1);
x2_cu  = x_cu(:,2)/scaleX_in_g(2,1);
x3_cu  = x_cu(:,3)/scaleX_in_g(3,1);
x4_cu  = x_cu(:,4)/scaleX_in_g(4,1);
x5_cu  = x_cu(:,5)/scaleX_in_g(5,1);
x6_cu  = x_cu(:,6)/scaleX_in_g(6,1);
 
Egfunc_cu = g0(1,1)+gx(1,1).*x1_cu+gx(1,2).*x2_cu+gx(1,3).*x3_cu+gx(1,4).*x4_cu+gx(1,5).*x5_cu+gx(1,6).*x6_cu;
if appOrder > 1
    Egfunc_cu = Egfunc_cu+gxxV(1,1).*x1_cu.*x1_cu+2.*gxxV(1,2).*x1_cu.*x2_cu+2.*gxxV(1,3).*x1_cu.*x3_cu+2.*gxxV(1,4).*x1_cu.*x4_cu+2.*gxxV(1,5).*x1_cu.*x5_cu+2.*gxxV(1,6).*x1_cu.*x6_cu+gxxV(1,7).*x2_cu.*x2_cu+2.*gxxV(1,8).*x2_cu.*x3_cu+2.*gxxV(1,9).*x2_cu.*x4_cu+2.*gxxV(1,10).*x2_cu.*x5_cu+2.*gxxV(1,11).*x2_cu.*x6_cu+gxxV(1,12).*x3_cu.*x3_cu+2.*gxxV(1,13).*x3_cu.*x4_cu+2.*gxxV(1,14).*x3_cu.*x5_cu+2.*gxxV(1,15).*x3_cu.*x6_cu+gxxV(1,16).*x4_cu.*x4_cu+2.*gxxV(1,17).*x4_cu.*x5_cu+2.*gxxV(1,18).*x4_cu.*x6_cu+gxxV(1,19).*x5_cu.*x5_cu+2.*gxxV(1,20).*x5_cu.*x6_cu+gxxV(1,21).*x6_cu.*x6_cu;
end
if appOrder > 2
    Egfunc_cu = Egfunc_cu+gxxxV(1,1).*x1_cu.*x1_cu.*x1_cu+3.*gxxxV(1,2).*x1_cu.*x1_cu.*x2_cu+3.*gxxxV(1,3).*x1_cu.*x1_cu.*x3_cu+3.*gxxxV(1,4).*x1_cu.*x1_cu.*x4_cu+3.*gxxxV(1,5).*x1_cu.*x1_cu.*x5_cu+3.*gxxxV(1,6).*x1_cu.*x1_cu.*x6_cu+3.*gxxxV(1,7).*x1_cu.*x2_cu.*x2_cu+6.*gxxxV(1,8).*x1_cu.*x2_cu.*x3_cu+6.*gxxxV(1,9).*x1_cu.*x2_cu.*x4_cu+6.*gxxxV(1,10).*x1_cu.*x2_cu.*x5_cu+6.*gxxxV(1,11).*x1_cu.*x2_cu.*x6_cu+3.*gxxxV(1,12).*x1_cu.*x3_cu.*x3_cu+6.*gxxxV(1,13).*x1_cu.*x3_cu.*x4_cu+6.*gxxxV(1,14).*x1_cu.*x3_cu.*x5_cu+6.*gxxxV(1,15).*x1_cu.*x3_cu.*x6_cu+3.*gxxxV(1,16).*x1_cu.*x4_cu.*x4_cu+6.*gxxxV(1,17).*x1_cu.*x4_cu.*x5_cu+6.*gxxxV(1,18).*x1_cu.*x4_cu.*x6_cu+3.*gxxxV(1,19).*x1_cu.*x5_cu.*x5_cu+6.*gxxxV(1,20).*x1_cu.*x5_cu.*x6_cu+3.*gxxxV(1,21).*x1_cu.*x6_cu.*x6_cu+gxxxV(1,22).*x2_cu.*x2_cu.*x2_cu+3.*gxxxV(1,23).*x2_cu.*x2_cu.*x3_cu+3.*gxxxV(1,24).*x2_cu.*x2_cu.*x4_cu+3.*gxxxV(1,25).*x2_cu.*x2_cu.*x5_cu+3.*gxxxV(1,26).*x2_cu.*x2_cu.*x6_cu+3.*gxxxV(1,27).*x2_cu.*x3_cu.*x3_cu+6.*gxxxV(1,28).*x2_cu.*x3_cu.*x4_cu+6.*gxxxV(1,29).*x2_cu.*x3_cu.*x5_cu+6.*gxxxV(1,30).*x2_cu.*x3_cu.*x6_cu+3.*gxxxV(1,31).*x2_cu.*x4_cu.*x4_cu+6.*gxxxV(1,32).*x2_cu.*x4_cu.*x5_cu+6.*gxxxV(1,33).*x2_cu.*x4_cu.*x6_cu+3.*gxxxV(1,34).*x2_cu.*x5_cu.*x5_cu+6.*gxxxV(1,35).*x2_cu.*x5_cu.*x6_cu+3.*gxxxV(1,36).*x2_cu.*x6_cu.*x6_cu+gxxxV(1,37).*x3_cu.*x3_cu.*x3_cu+3.*gxxxV(1,38).*x3_cu.*x3_cu.*x4_cu+3.*gxxxV(1,39).*x3_cu.*x3_cu.*x5_cu+3.*gxxxV(1,40).*x3_cu.*x3_cu.*x6_cu+3.*gxxxV(1,41).*x3_cu.*x4_cu.*x4_cu+6.*gxxxV(1,42).*x3_cu.*x4_cu.*x5_cu+6.*gxxxV(1,43).*x3_cu.*x4_cu.*x6_cu+3.*gxxxV(1,44).*x3_cu.*x5_cu.*x5_cu+6.*gxxxV(1,45).*x3_cu.*x5_cu.*x6_cu+3.*gxxxV(1,46).*x3_cu.*x6_cu.*x6_cu+gxxxV(1,47).*x4_cu.*x4_cu.*x4_cu+3.*gxxxV(1,48).*x4_cu.*x4_cu.*x5_cu+3.*gxxxV(1,49).*x4_cu.*x4_cu.*x6_cu+3.*gxxxV(1,50).*x4_cu.*x5_cu.*x5_cu+6.*gxxxV(1,51).*x4_cu.*x5_cu.*x6_cu+3.*gxxxV(1,52).*x4_cu.*x6_cu.*x6_cu+gxxxV(1,53).*x5_cu.*x5_cu.*x5_cu+3.*gxxxV(1,54).*x5_cu.*x5_cu.*x6_cu+3.*gxxxV(1,55).*x5_cu.*x6_cu.*x6_cu+gxxxV(1,56).*x6_cu.*x6_cu.*x6_cu;
end
 
%% The control variable in the next period 
xVec_cup = zeros(numX,nx);
mx = setup.mx;
if mx > 0 
end 
if setup.myx > 0 
    xVec_cup(:,mx+1) = c_cu(:,1)-gSS(1,1);
end 
diaghx = diag(setup.hx);
xVec_cup(:,nx1+1:nx)  = x_cu(:,nx1+1:nx).*repmat(transpose(diaghx(nx1+1:nx)),numX,1); 
  
%Adding shocks to the states at t + 1
diagEta = diag(setup.eta(nx1+1:nx,:));
x1_cup  = repmat(xVec_cup(:,1),1,nNodes);
x2_cup  = repmat(xVec_cup(:,2),1,nNodes) +2./(exp(-OMEGAz.*x_cu(:,2))+1).*diagEta(1,1).*repmat(epsNodes(1,:),numX,1);
x3_cup  = repmat(xVec_cup(:,3),1,nNodes) +2./(exp(-OMEGAd.*x_cu(:,3))+1).*diagEta(2,1).*repmat(epsNodes(2,:),numX,1);
x4_cup  = repmat(xVec_cup(:,4),1,nNodes) +2./(exp(-OMEGAn.*x_cu(:,4))+1).*diagEta(3,1).*repmat(epsNodes(3,:),numX,1);
x5_cup  = repmat(xVec_cup(:,5),1,nNodes) +2./(exp(-OMEGAp.*x_cu(:,5))+1).*diagEta(4,1).*repmat(epsNodes(4,:),numX,1);
x6_cup  = repmat(xVec_cup(:,6),1,nNodes) +2./(exp(-OMEGAa.*x_cu(:,6))+1).*diagEta(5,1).*repmat(epsNodes(5,:),numX,1);
 
% Scaling the states at time t+1 by the scaleX_in_g
x1_cup  = x1_cup/scaleX_in_g(1,1);
x2_cup  = x2_cup/scaleX_in_g(2,1);
x3_cup  = x3_cup/scaleX_in_g(3,1);
x4_cup  = x4_cup/scaleX_in_g(4,1);
x5_cup  = x5_cup/scaleX_in_g(5,1);
x6_cup  = x6_cup/scaleX_in_g(6,1);
 
var_cup = var0(1,1)+varx(1,1).*x1_cup+varx(1,2).*x2_cup+varx(1,3).*x3_cup+varx(1,4).*x4_cup+varx(1,5).*x5_cup+varx(1,6).*x6_cup;
if appOrder > 1
    var_cup = var_cup+varxx(1,1).*x1_cup.*x1_cup+2.*varxx(1,2).*x1_cup.*x2_cup+2.*varxx(1,3).*x1_cup.*x3_cup+2.*varxx(1,4).*x1_cup.*x4_cup+2.*varxx(1,5).*x1_cup.*x5_cup+2.*varxx(1,6).*x1_cup.*x6_cup+varxx(1,7).*x2_cup.*x2_cup+2.*varxx(1,8).*x2_cup.*x3_cup+2.*varxx(1,9).*x2_cup.*x4_cup+2.*varxx(1,10).*x2_cup.*x5_cup+2.*varxx(1,11).*x2_cup.*x6_cup+varxx(1,12).*x3_cup.*x3_cup+2.*varxx(1,13).*x3_cup.*x4_cup+2.*varxx(1,14).*x3_cup.*x5_cup+2.*varxx(1,15).*x3_cup.*x6_cup+varxx(1,16).*x4_cup.*x4_cup+2.*varxx(1,17).*x4_cup.*x5_cup+2.*varxx(1,18).*x4_cup.*x6_cup+varxx(1,19).*x5_cup.*x5_cup+2.*varxx(1,20).*x5_cup.*x6_cup+varxx(1,21).*x6_cup.*x6_cup;
end
if appOrder > 2
    var_cup = var_cup+varxxx(1,1).*x1_cup.*x1_cup.*x1_cup+3.*varxxx(1,2).*x1_cup.*x1_cup.*x2_cup+3.*varxxx(1,3).*x1_cup.*x1_cup.*x3_cup+3.*varxxx(1,4).*x1_cup.*x1_cup.*x4_cup+3.*varxxx(1,5).*x1_cup.*x1_cup.*x5_cup+3.*varxxx(1,6).*x1_cup.*x1_cup.*x6_cup+3.*varxxx(1,7).*x1_cup.*x2_cup.*x2_cup+6.*varxxx(1,8).*x1_cup.*x2_cup.*x3_cup+6.*varxxx(1,9).*x1_cup.*x2_cup.*x4_cup+6.*varxxx(1,10).*x1_cup.*x2_cup.*x5_cup+6.*varxxx(1,11).*x1_cup.*x2_cup.*x6_cup+3.*varxxx(1,12).*x1_cup.*x3_cup.*x3_cup+6.*varxxx(1,13).*x1_cup.*x3_cup.*x4_cup+6.*varxxx(1,14).*x1_cup.*x3_cup.*x5_cup+6.*varxxx(1,15).*x1_cup.*x3_cup.*x6_cup+3.*varxxx(1,16).*x1_cup.*x4_cup.*x4_cup+6.*varxxx(1,17).*x1_cup.*x4_cup.*x5_cup+6.*varxxx(1,18).*x1_cup.*x4_cup.*x6_cup+3.*varxxx(1,19).*x1_cup.*x5_cup.*x5_cup+6.*varxxx(1,20).*x1_cup.*x5_cup.*x6_cup+3.*varxxx(1,21).*x1_cup.*x6_cup.*x6_cup+varxxx(1,22).*x2_cup.*x2_cup.*x2_cup+3.*varxxx(1,23).*x2_cup.*x2_cup.*x3_cup+3.*varxxx(1,24).*x2_cup.*x2_cup.*x4_cup+3.*varxxx(1,25).*x2_cup.*x2_cup.*x5_cup+3.*varxxx(1,26).*x2_cup.*x2_cup.*x6_cup+3.*varxxx(1,27).*x2_cup.*x3_cup.*x3_cup+6.*varxxx(1,28).*x2_cup.*x3_cup.*x4_cup+6.*varxxx(1,29).*x2_cup.*x3_cup.*x5_cup+6.*varxxx(1,30).*x2_cup.*x3_cup.*x6_cup+3.*varxxx(1,31).*x2_cup.*x4_cup.*x4_cup+6.*varxxx(1,32).*x2_cup.*x4_cup.*x5_cup+6.*varxxx(1,33).*x2_cup.*x4_cup.*x6_cup+3.*varxxx(1,34).*x2_cup.*x5_cup.*x5_cup+6.*varxxx(1,35).*x2_cup.*x5_cup.*x6_cup+3.*varxxx(1,36).*x2_cup.*x6_cup.*x6_cup+varxxx(1,37).*x3_cup.*x3_cup.*x3_cup+3.*varxxx(1,38).*x3_cup.*x3_cup.*x4_cup+3.*varxxx(1,39).*x3_cup.*x3_cup.*x5_cup+3.*varxxx(1,40).*x3_cup.*x3_cup.*x6_cup+3.*varxxx(1,41).*x3_cup.*x4_cup.*x4_cup+6.*varxxx(1,42).*x3_cup.*x4_cup.*x5_cup+6.*varxxx(1,43).*x3_cup.*x4_cup.*x6_cup+3.*varxxx(1,44).*x3_cup.*x5_cup.*x5_cup+6.*varxxx(1,45).*x3_cup.*x5_cup.*x6_cup+3.*varxxx(1,46).*x3_cup.*x6_cup.*x6_cup+varxxx(1,47).*x4_cup.*x4_cup.*x4_cup+3.*varxxx(1,48).*x4_cup.*x4_cup.*x5_cup+3.*varxxx(1,49).*x4_cup.*x4_cup.*x6_cup+3.*varxxx(1,50).*x4_cup.*x5_cup.*x5_cup+6.*varxxx(1,51).*x4_cup.*x5_cup.*x6_cup+3.*varxxx(1,52).*x4_cup.*x6_cup.*x6_cup+varxxx(1,53).*x5_cup.*x5_cup.*x5_cup+3.*varxxx(1,54).*x5_cup.*x5_cup.*x6_cup+3.*varxxx(1,55).*x5_cup.*x6_cup.*x6_cup+varxxx(1,56).*x6_cup.*x6_cup.*x6_cup;
end
 
% Numerical integration
Evar_cu = var_cup*wNodes;
% The Euler errors 
resEuler = -Egfunc_cu + Evar_cu;
resEulerWeight = resEuler.*(1./(1+sum(x_cu.^2,2)));
res = resEulerWeight(:);
end
